#include "sgPetscSpace.h"
// The PETSC_FUNCTION_NAME_CXX will be defined only if the '--with-cxx' is enabled.
// However, it seems that the PETSc installed by pacman was built without this
// configuration.
// see https://petsc.org/release/install/install/
// see petscconfig.h

// #define PETSC_FUNCTION_NAME_CXX __func__
// #define PETSC_FUNCTION_NAME_C __func__

// #include "petscksp.h"
// #include "petscpc.h"
// #include "petscsys.h"
#include <iostream>

#include <petsc.h>
// #include "omp.h"
#include "WrappedPetsc.h"

using namespace SG::Algebra;

namespace  // Anonymous
{
    // Convertion between PetscScalar and builtin types.
    template <typename TDataType>
    TDataType P2C (const PetscScalar& petsc_value);

    template <>
    Real_t P2C<Real_t> (const PetscScalar& petsc_value)
    {
#if defined(PETSC_HAVE_COMPLEX)
        return PetscRealPart (petsc_value);
#else
        return petsc_value;
#endif
    }

    template <>
    Complex_t P2C<Complex_t> (const PetscScalar& petsc_value)
    {
#if defined(PETSC_HAVE_COMPLEX)
        return Complex_t (PetscRealPart (petsc_value), PetscImaginaryPart (petsc_value));
#else
        return Complex_t (petsc_value, 0.0);
#endif
    }

    template <typename TDataType>
    PetscScalar C2P (const TDataType& petsc_value);

    template <>
    PetscScalar C2P<Real_t> (const Real_t& value)
    {
        PetscScalar petsc_value;
#if defined(PETSC_HAVE_COMPLEX)
        petsc_value = PetscCMPLX (value, 0.0);
        return petsc_value;
#else
        petsc_value = value;
        return petsc_value;
#endif
    }

    template <>
    PetscScalar C2P<Complex_t> (const Complex_t& value)
    {
        PetscScalar petsc_value;
#if defined(PETSC_HAVE_COMPLEX)
        petsc_value = PetscCMPLX (value.real (), value.imag ());
        return petsc_value;
#else
        petsc_value = value.real ();
        return petsc_value;
#endif
    }

    // /// Copy blocks from Vec into local vectors
    // template <typename TDataType>
    // std::vector<TDataType> get_vector_values (const Vec x, const IndexArray_t& indices)
    // {
    //     std::vector<TDataType> c_values;
    //     PetscInt               n     = indices.size ();

    //     // src vector and index set
    //     Vec      src           = x;
    //     std::vector<PetscInt>  src_indices;
    //     for (auto i = 0; i < indices.size (); ++ i)
    //     {
    //         src_indices.push_back (indices[i]);
    //     }

    //     MPI_Comm src_comm;
    //     PetscObjectGetComm (( PetscObject ) x, &src_comm);
    //     IS src_is;  // if NULL scatters all values
    //     ISCreateGeneral (src_comm, n, src_indices.data (), PETSC_USE_POINTER, &src_is);

    //     // dest vector and index set
    //     Vec dest;
    //     VecCreateSeq (PETSC_COMM_SELF, 1, &dest);
    //     IS dest_is = PETSC_NULLPTR;  // if NULL fills entire vector yin in order

    //     // vector scatter
    //     VecScatter scatter;
    //     VecScatterCreate (src, /*src is=*/src_is, /*dest vec=*/dest, /*dest is=*/dest_is, &scatter);
    //     // Actually perform the scatter
    //     VecScatterBegin (scatter, src, dest, INSERT_VALUES, SCATTER_FORWARD);
    //     VecScatterEnd (scatter, src, dest, INSERT_VALUES, SCATTER_FORWARD);

    //     PetscScalar* values;
    //     VecGetArray (dest, &values);
    //     for (auto i = 0; i < n; ++ i)
    //     {
    //         PetscScalar petsc_value = values[0];
    //         TDataType   value       = P2C<TDataType> (petsc_value);
    //         c_values.push_back (value);
    //     }
    //     VecRestoreArray (dest, &values);

    //     // VecDestroy (&src);
    //     ISDestroy (&src_is);
    //     VecDestroy (&dest);
    //     ISDestroy (&dest_is);
    //     VecScatterDestroy (&scatter);

    //     return c_values;
    // }

    /// Get the specified off-process value from the parallel vector.
    template <typename TDataType>
    TDataType get_vector_value (const Vec x, const Index_t& index)
    {
        // src vector and index set
        Vec      src           = x;
        PetscInt src_indices[] = { index };

        MPI_Comm src_comm;
        PetscObjectGetComm (( PetscObject ) x, &src_comm);
        IS src_is;  // if NULL scatters all values
        ISCreateGeneral (src_comm, 1, src_indices, PETSC_USE_POINTER, &src_is);

        // dest vector and index set
        Vec dest;
        VecCreateSeq (PETSC_COMM_SELF, 1, &dest);

        IS dest_is = PETSC_NULLPTR;  // if NULL fills entire vector yin in order

        // vector scatter
        VecScatter scatter;
        VecScatterCreate (src, /*src is=*/src_is, /*dest vec=*/dest, /*dest is=*/dest_is, &scatter);
        // Actually perform the scatter
        VecScatterBegin (scatter, src, dest, INSERT_VALUES, SCATTER_FORWARD);
        VecScatterEnd (scatter, src, dest, INSERT_VALUES, SCATTER_FORWARD);

        PetscScalar* values;
        VecGetArray (dest, &values);
        PetscScalar petsc_value = values[0];
        TDataType   value       = P2C<TDataType> (petsc_value);
        VecRestoreArray (dest, &values);

        // VecDestroy (&src);
        ISDestroy (&src_is);
        VecDestroy (&dest);
        ISDestroy (&dest_is);
        VecScatterDestroy (&scatter);

        return value;
    }

    /// Get the specified off-process value from the parallel matrix.
    template <typename TDataType>
    TDataType get_matrix_value (const Mat mat, const Index_t& i, const Index_t& j)
    {
        // Gets the values from a given column of a matrix.
        // The vector must have the same parallel row layout as the matrix.
        MPI_Comm comm;
        PetscObjectGetComm (( PetscObject ) mat, &comm);

        Vec      yy;
        PetscInt col = j;
        VecCreate (comm, &yy);
        CHKERRQ (MatGetColumnVector (mat, yy, col));
        auto value = get_vector_value<TDataType> (yy, i);

        return value;
    }
}

namespace SG
{
    namespace Algebra
    {
        struct PetscVectorPrivate
        {
            Vec      vec = nullptr;
            PetscInt n   = PETSC_DECIDE;     // the local size (or PETSC_DECIDE to have it set)
            PetscInt N   = PETSC_DETERMINE;  // the global size (or PETSC_DETERMINE to have it set)
        };

        struct PetscDenseMatrixPrivate
        {
            Mat mat = nullptr;
        };

        struct PetscSparseMatrixPrivate
        {
            Mat      mat = nullptr;
            PetscInt m   = PETSC_DECIDE;  // number of local rows (or PETSC_DECIDE)
            PetscInt n   = PETSC_DECIDE;  // number of local columns (or PETSC_DECIDE)
            PetscInt M   = PETSC_DETERMINE;  // number of global rows (or PETSC_DETERMINE)
            PetscInt N   = PETSC_DETERMINE;  // number of global columns (or PETSC_DETERMINE)
        };

    }  // namespace Algebra
}  // namespace SG

//=============================================================================
// PetscVector::Entry

template <typename TDataType>
typename PetscVector<TDataType>::Entry& PetscVector<TDataType>::Entry::operator= (const PetscVector<TDataType>::Entry& r)
{
    // as explained in the class
    // documentation, this is not the copy
    // operator. so simply pass on to the
    // "correct" assignment operator
    *this = static_cast<TDataType> (r);

    return *this;
}

template <typename TDataType>
const typename PetscVector<TDataType>::Entry& PetscVector<TDataType>::Entry::operator= (const TDataType& value) const
{
    const PetscInt    petsc_i     = mIndex;
    const PetscScalar petsc_value = C2P<TDataType> (value);
    VecSetValues (mVector, 1, &petsc_i, &petsc_value, INSERT_VALUES);

    return *this;
}

template <typename TDataType>
const typename PetscVector<TDataType>::Entry& PetscVector<TDataType>::Entry::operator+= (const TDataType& value) const
{
    // use the PETSc function to add something
    const PetscInt    petsc_i     = mIndex;
    const PetscScalar petsc_value = C2P<TDataType> (value);
    VecSetValues (mVector, 1, &petsc_i, &petsc_value, ADD_VALUES);

    return *this;
}

template <typename TDataType>
const typename PetscVector<TDataType>::Entry& PetscVector<TDataType>::Entry::operator-= (const TDataType& value) const
{
    // use the PETSc function to
    // add something
    const PetscInt    petsc_i     = mIndex;
    const PetscScalar petsc_value = C2P<TDataType> (-value);
    VecSetValues (mVector, 1, &petsc_i, &petsc_value, ADD_VALUES);

    return *this;
}

template <typename TDataType>
const typename PetscVector<TDataType>::Entry& PetscVector<TDataType>::Entry::operator*= (const TDataType& value) const
{
    const PetscInt    petsc_i   = mIndex;
    const PetscScalar new_value = C2P<TDataType> (static_cast<TDataType> (*this) * value);
    VecSetValues (mVector, 1, &petsc_i, &new_value, INSERT_VALUES);

    return *this;
}

template <typename TDataType>
const typename PetscVector<TDataType>::Entry& PetscVector<TDataType>::Entry::operator/= (const TDataType& value) const
{
    const PetscInt    petsc_i   = mIndex;
    const PetscScalar new_value = static_cast<TDataType> (*this) / value;
    VecSetValues (mVector, 1, &petsc_i, &new_value, INSERT_VALUES);

    return *this;
}

template <typename TDataType>
PetscVector<TDataType>::Entry::operator TDataType () const
{
    // first verify that the requested
    // element is actually locally
    // available
    const PetscInt petsc_i     = mIndex;
    const PetscInt petsc_begin = mVector.row_start ();
    const PetscInt petsc_end   = mVector.row_end ();

    if (petsc_i >= petsc_begin && petsc_i < petsc_end)
    {
        PetscScalar petsc_value;
        VecGetValues (( Vec ) mVector, 1, &petsc_i, &petsc_value);
        TDataType value = P2C<TDataType> (petsc_value);

        return value;
    }

    /*
     * VecGetValues() can only access local values from the vector. To get of-process values,
     * the user should create a new vector where the components will be stored and then perform
     *  the appropriate vector scatter.
     */

    // assert (0
    //         && "You are trying to access an element of a vector "
    //            "that is neither a locally owned element nor a "
    //            "ghost element of the vector.");

#ifdef SG_VERBOSE
    std::cerr << "You are trying to access an element of a vector "
              << "that is neither a locally owned element nor a "
              << "ghost element of the vector." << std::endl;
#endif

    Vec  src   = ( Vec ) mVector;
    auto value = get_vector_value<TDataType> (src, mIndex);

    return value;
}

template <typename TDataType>
Real_t PetscVector<TDataType>::Entry::real () const
{
    TDataType value = static_cast<TDataType> (*this);

    return std::real (value);
}

template <typename TDataType>
Real_t PetscVector<TDataType>::Entry::imag () const
{
    TDataType value = static_cast<TDataType> (*this);

    return std::imag (value);
}

//=============================================================================
// PetscVector

template <typename TDataType>
PetscVector<TDataType>::PetscVector ()
{
    d = new PetscVectorPrivate ();
    d->vec = PETSC_NULLPTR;

    // CHKERRV (VecCreate (PETSC_COMM_WORLD, &d->vec));
    // this->mInitialized = true;
}

template <typename TDataType>
PetscVector<TDataType>::PetscVector (Sparsity* sp)
{
    d = new PetscVectorPrivate ();

    // Because PETSC types are generally distributed across, and accessible from,
    // multiple MPI pro-cesses, the ﬁrst argument of an ObjectCreate() method is
    // an MPI communicator (“COMM”).
    // see PETSc for Partial Differential Equations
    CHKERRV (VecCreate (PETSC_COMM_WORLD, &d->vec));
    // The type of storage for the vector may be set with either calls to VecSetType()
    // or VecSetFromOptions()
    CHKERRV (VecSetUp (d->vec));
    // CHKERRV (VecSetSizes (d->vec, d->n, d->N));

    // this->mInitialized = true;
}

template <typename TDataType>
PetscVector<TDataType>::PetscVector (const Index_t& N)
{
    d = new PetscVectorPrivate ();
    d->n = PETSC_DECIDE;
    d->N = N;

    // Because PETSC types are generally distributed across, and accessible from,
    // multiple MPI pro-cesses, the ﬁrst argument of an ObjectCreate() method is
    // an MPI communicator (“COMM”).
    // see PETSc for Partial Differential Equations
    CHKERRV (VecCreate (PETSC_COMM_WORLD, &d->vec));
    // The type of storage for the vector may be set with either calls to VecSetType()
    // or VecSetFromOptions()
    CHKERRV (VecSetSizes (d->vec, PETSC_DECIDE, N));
    CHKERRV (VecSetUp (d->vec));

    // this->mInitialized = true;
}

template <typename TDataType>
PetscVector<TDataType>::PetscVector (Vec vec, bool inc_ref_count)
{
    d      = new PetscVectorPrivate ();

    d->vec = vec;
    if (inc_ref_count)
    {
        PetscObjectReference (( PetscObject ) d->vec);
    }
}

template <typename TDataType>
PetscVector<TDataType>::PetscVector (const PetscVector<TDataType>& x)
{
    d = new PetscVectorPrivate ();

    CHKERRV (VecDuplicate (x.d->vec, &d->vec));
    CHKERRV (VecCopy (x.d->vec, d->vec));

    // this->mInitialized = true;
}

template <typename TDataType>
PetscVector<TDataType>::PetscVector (PetscVector<TDataType>&& v)
{
    d      = new PetscVectorPrivate ();
    d->vec = std::exchange (v.d->vec, nullptr);

    // this->mInitialized = true;
}

template <typename TDataType>
PetscVector<TDataType>::~PetscVector ()
{
    if (d && d->vec)
    {
        // If we encounter an error here, print a warning but otherwise
        // keep going since we may be recovering from an exception.
        CHKERRV (VecDestroy (&d->vec));
        d->vec = nullptr;
        // this->mInitialized = false;
    }

    if (d)
    {
        delete d;
    }
}

template <typename TDataType>
PetscVector<TDataType>& PetscVector<TDataType>::operator= (const PetscVector<TDataType>& v)
{
    VecDuplicate (v.d->vec, &d->vec);
    VecCopy (v.d->vec, d->vec);

    return *this;
}

template <typename TDataType>
PetscVector<TDataType>& PetscVector<TDataType>::operator= (PetscVector<TDataType>&& v)
{
    std::swap (d->vec, v.d->vec);
    return *this;
}

template <typename TDataType>
PetscVector<TDataType>::operator const Vec& () const
{
    return d->vec;
}

template <typename TDataType>
PetscVector<TDataType>::operator Vec& ()
{
    return d->vec;
}

template <typename TDataType>
Vec PetscVector<TDataType>::vec ()
{
    return d->vec;
}

template <typename TDataType>
Vec PetscVector<TDataType>::vec () const
{
    return d->vec;
}

template <typename TDataType>
PetscVector<TDataType> PetscVector<TDataType>::clone () const
{
    Vec vec;
    VecDuplicate (d->vec, &vec);
    VecCopy (d->vec, vec);
    PetscVector<TDataType> y (vec, true);
    VecDestroy (&vec);
    return y;
}

template <typename TDataType>
Space<TDataType>* PetscVector<TDataType>::getSpace ()
{
    static PetscSpace<TDataType> petsc_space;
    return &petsc_space;
}

// template <typename TDataType>
// void PetscVector<TDataType>::init (Sparsity* sp)
// {
//     // Because PETSC types are generally distributed across, and accessible from,
//     // multiple MPI pro-cesses, the ﬁrst argument of an ObjectCreate() method is
//     // an MPI communicator (“COMM”).
//     // see PETSc for Partial Differential Equations
//     CHKERRV (VecCreate (PETSC_COMM_WORLD, &d->vec));
//     // CHKERRV (VecSetType (d->vec, VECMPI));
//     CHKERRV (VecSetSizes (d->vec, d->n, d->N));

//     this->mInitialized = true;
// }

template <typename TDataType>
void PetscVector<TDataType>::close ()
{
    VecAssemblyBegin (d->vec);
    VecAssemblyEnd (d->vec);
}

template <typename TDataType>
void PetscVector<TDataType>::clear ()
{
    if (d && d->vec)
    {
        // If we encounter an error here, print a warning but otherwise
        // keep going since we may be recovering from an exception.
        CHKERRV (VecDestroy (&d->vec));
        d->vec = nullptr;
        // this->mInitialized = false;
    }
}

template <typename TDataType>
void PetscVector<TDataType>::beginAssembly ()
{
    VecAssemblyBegin (d->vec);
}

template <typename TDataType>
void PetscVector<TDataType>::endAssembly ()
{
    VecAssemblyEnd (d->vec);
}

template <typename TDataType>
void PetscVector<TDataType>::assemble ()
{
    VecAssemblyBegin (d->vec);
    VecAssemblyEnd (d->vec);
}

// template <typename TDataType>
// bool PetscVector<TDataType>::assembled ()
// {
//     PetscBool assembled;
//     CHKERRQ (VecAssembled (d->vec, &assembled));

//     return assembled == PETSC_TRUE ? true : false;
// }

template <typename TDataType>
void PetscVector<TDataType>::resize (Index_t N)
{
    // has the exact dimensions?
    if (d->vec != PETSC_NULLPTR)
    {
        PetscInt N1;
        CHKERRV (VecGetSize (d->vec, &N1));

        if (N1 == N)
        {
            return;
        }
    }

    if (d->vec != PETSC_NULLPTR)
    {
        CHKERRV (VecDestroy (&d->vec));
    }

    CHKERRV (VecCreate (PETSC_COMM_WORLD, &d->vec));
    CHKERRV (VecSetSizes (d->vec, PETSC_DECIDE, N));
    CHKERRV (VecSetUp (d->vec));
}

template <typename TDataType>
void PetscVector<TDataType>::conservativeResize (Index_t N)
{
    PetscInt new_size = N;
    if (d->vec == PETSC_NULLPTR)
    {
        // d->N = N;
        CHKERRV (VecCreate (PETSC_COMM_WORLD, &d->vec));
        CHKERRV (VecSetSizes (d->vec, PETSC_DECIDE, N));
        CHKERRV (VecSetUp (d->vec));

        return;
    }

    // once the vector or matrices sizes have been set and the matrices or vectors
    // are fully usable one cannot change the size of the matrices or vectors or
    // number of processors they live on. One may create new vectors and copy,
    // for example using VecScatterCreate(), the old values from the previous vector.

    MPI_Comm comm;
    PetscObjectGetComm (( PetscObject ) d->vec, &comm);
    IS  src_is = PETSC_NULLPTR;  // if NULL scatters all values
    Vec src    = d->vec;

    PetscInt old_size;
    CHKERRV (VecGetSize (d->vec, &old_size));
    IS dest_is = PETSC_NULLPTR;  // if NULL fills entire vector yin in order
    CHKERRV (ISCreateStride (comm, old_size, 0, 1, &dest_is));

    // dest vector and index set
    Vec dest;
    CHKERRV (VecCreate (comm, &dest));
    CHKERRV (VecSetSizes (dest, PETSC_DECIDE, new_size));
    CHKERRV (VecSetUp (dest));

    // vector scatter
    VecScatter scatter;
    CHKERRV (VecScatterCreate (src, src_is, dest, dest_is, &scatter));
    // Actually perform the scatter
    CHKERRV (VecScatterBegin (scatter, src, dest, INSERT_VALUES, SCATTER_FORWARD));
    CHKERRV (VecScatterEnd (scatter, src, dest, INSERT_VALUES, SCATTER_FORWARD));
    CHKERRV (VecDestroy (&d->vec));
    d->vec = dest;
}

template <typename TDataType>
Index_t PetscVector<TDataType>::rows () const
{
    PetscInt petsc_size = 0;

    CHKERRQ (VecGetSize (d->vec, &petsc_size));

    return static_cast<Index_t> (petsc_size);
}

template <typename TDataType>
Index_t PetscVector<TDataType>::row_start () const
{
    Index_t first = 0;

    PetscInt petsc_first = 0, petsc_last = 0;
    CHKERRQ (VecGetOwnershipRange (d->vec, &petsc_first, &petsc_last));
    first = static_cast<Index_t> (petsc_first);

    return first;
}

template <typename TDataType>
Index_t PetscVector<TDataType>::row_end () const
{
    Index_t last = 0;

    PetscInt petsc_first = 0, petsc_last = 0;
    CHKERRQ (VecGetOwnershipRange (d->vec, &petsc_first, &petsc_last));
    last = static_cast<Index_t> (petsc_last);

    return last;
}

template <typename TDataType>
void PetscVector<TDataType>::setZero ()
{
    CHKERRV (VecZeroEntries (d->vec));
}

template <typename TDataType>
void PetscVector<TDataType>::set (const Index_t& i, const TDataType& value)
{
    PetscInt    i_val       = static_cast<PetscInt> (i);
    PetscScalar petsc_value = value;

    CHKERRV (VecSetValues (d->vec, 1, &i_val, &petsc_value, INSERT_VALUES));
}

template <typename TDataType>
void PetscVector<TDataType>::add (const Index_t& i, const TDataType& value)
{
    PetscInt    i_val       = i;
    PetscScalar petsc_value = value;

    CHKERRV (VecSetValues (d->vec, 1, &i_val, &petsc_value, ADD_VALUES));
}

template <typename TDataType>
typename PetscVector<TDataType>::Entry PetscVector<TDataType>::operator() (const Index_t& i)
{
    Entry e (*this, i);
    return e;
}

template <typename TDataType>
TDataType PetscVector<TDataType>::operator() (const Index_t& i) const
{
    // first verify that the requested
    // element is actually locally
    // available
    PetscInt petsc_i = i;
    PetscInt petsc_begin, petsc_end;

    CHKERRQ (VecGetOwnershipRange (d->vec, &petsc_begin, &petsc_end));

    if (petsc_i >= petsc_begin && petsc_i < petsc_end)
    {
        PetscScalar petsc_value;
        VecGetValues (d->vec, 1, &petsc_i, &petsc_value);
        TDataType value = P2C<TDataType> (petsc_value);

        return value;
    }

    /*
     * VecGetValues() can only access local values from the vector. To get of-process values,
     * the user should create a new vector where the components will be stored and then perform
     *  the appropriate vector scatter.
     */

    // assert (0
    //         && "You are trying to access an element of a vector "
    //            "that is neither a locally owned element nor a "
    //            "ghost element of the vector.");

#ifdef SG_VERBOSE
    std::cerr << "You are trying to access an element of a vector "
              << "that is neither a locally owned element nor a "
              << "ghost element of the vector." << std::endl;
#endif

    Vec  src   = d->vec;
    auto value = get_vector_value<TDataType> (src, i);

    return value;
}

template <typename TDataType>
typename PetscVector<TDataType>::Entry PetscVector<TDataType>::operator[] (const Index_t& i)
{
    return this->operator()(i);
}

template <typename TDataType>
TDataType PetscVector<TDataType>::operator[] (const Index_t& i) const
{
    return this->operator()(i);
}

template <typename TDataType>
void PetscVector<TDataType>::print (std::ostream& out) const
{
    /**
     * To print the vector to the screen, one can use the viewer PETSC_VIEWER_STDOUT_WORLD,
     * which ensures that parallel vectors are printed correctly to stdout. To display the
     * vector in an X-window, one can use the default X-windows viewer PETSC_VIEWER_DRAW_WORLD,
     * or one can create a viewer with the routine PetscViewerDrawOpen()
     */
    VecView (d->vec, PETSC_VIEWER_STDOUT_WORLD);
}

template <typename TDataType>
void PetscVector<TDataType>::create_subvector (Vector<TDataType>& subvector, const GeneralIndexSet& rows) const
{
    // Make sure the passed in subvector is really a PetscVector
    PetscVector<TDataType>* petsc_subvector = dynamic_cast<PetscVector<TDataType>*> (&subvector);
    petsc_subvector->clear ();

    MPI_Comm comm;
    CHKERRV (PetscObjectGetComm (( PetscObject ) d->vec, &comm));

    std::vector<PetscInt> idx;
    for (int i = 0; i < rows.size (); ++ i)
    {
        idx.push_back (rows[i]);
    }

    IS parent_is;
    CHKERRV (ISCreateGeneral (comm, idx.size (), idx.data (), PETSC_USE_POINTER, &parent_is));

    Vec subvec;
    CHKERRV (VecGetSubVector (d->vec, parent_is, &subvec));
    petsc_subvector->d->vec = subvec;

    ISDestroy (&parent_is);
}

template <typename TDataType>
void PetscVector<TDataType>::create_subvector (Vector<TDataType>& subvector, const StrideIndexSet& stride) const
{
    GeneralIndexSet idx;

    for (auto i = 0; i < stride.n; ++i)
    {
        idx.push_back (stride.first + i * stride.step);
    }

    create_subvector (subvector, idx);
}

template <typename TDataType>
void PetscVector<TDataType>::max (Index_t& index, TDataType& value) const
{
    PetscInt  petsc_index;
    PetscReal petsc_value;

    CHKERRV (VecMax (d->vec, &petsc_index, &petsc_value));

    index = petsc_index;
    value = petsc_value;
}

template <typename TDataType>
void PetscVector<TDataType>::linfty_norm (TDataType& value) const
{
    PetscReal petsc_value;
    CHKERRV (VecNorm (d->vec, NORM_INFINITY, &petsc_value));
    value = petsc_value;
}

//=============================================================================
// PetscDenseMatrix

template <typename TDataType>
PetscDenseMatrix<TDataType>::PetscDenseMatrix ()
{
    d = new PetscDenseMatrixPrivate ();
}

template <typename TDataType>
PetscDenseMatrix<TDataType>::PetscDenseMatrix (Index_t rows, Index_t cols)
{
    d = new PetscDenseMatrixPrivate ();
}

template <typename TDataType>
PetscDenseMatrix<TDataType>::~PetscDenseMatrix ()
{
    if (d)
    {
        delete d;
    }
}

template <typename TDataType>
void PetscDenseMatrix<TDataType>::print (std::ostream& out) const
{
    // CHKERRV (MatView (d->mat, PETSC_VIEWER_STDOUT_WORLD));
}

//=============================================================================
// PetscSparseMatrix::Entry

template <typename TDataType>
typename PetscSparseMatrix<TDataType>::Entry& PetscSparseMatrix<TDataType>::Entry::operator= (const PetscSparseMatrix<TDataType>::Entry& r)
{
    // as explained in the class
    // documentation, this is not the copy
    // operator. so simply pass on to the
    // "correct" assignment operator
    *this = static_cast<TDataType> (r);

    return *this;
}

template <typename TDataType>
const typename PetscSparseMatrix<TDataType>::Entry& PetscSparseMatrix<TDataType>::Entry::operator= (const TDataType& value) const
{
    const PetscInt    petsc_i     = mRow;
    const PetscInt    petsc_j     = mCol;
    const PetscScalar petsc_value = C2P<TDataType> (value);
    MatSetValues (( Mat ) mMatrix, 1, &petsc_i, 1, &petsc_j, &petsc_value, INSERT_VALUES);

    return *this;
}

template <typename TDataType>
const typename PetscSparseMatrix<TDataType>::Entry& PetscSparseMatrix<TDataType>::Entry::operator+= (const TDataType& value) const
{
    // use the PETSc function to add something
    const PetscInt    petsc_i     = mRow;
    const PetscInt    petsc_j     = mCol;
    const PetscScalar petsc_value = C2P<TDataType> (value);
    MatSetValues (( Mat ) mMatrix, 1, &petsc_i, 1, &petsc_j, &petsc_value, ADD_VALUES);

    return *this;
}

template <typename TDataType>
const typename PetscSparseMatrix<TDataType>::Entry& PetscSparseMatrix<TDataType>::Entry::operator-= (const TDataType& value) const
{
    // use the PETSc function to
    // add something
    const PetscInt    petsc_i     = mRow;
    const PetscInt    petsc_j     = mCol;
    const PetscScalar petsc_value = C2P<TDataType> (-value);
    MatSetValues (( Mat ) mMatrix, 1, &petsc_i, 1, &petsc_j, &petsc_value, ADD_VALUES);

    return *this;
}

template <typename TDataType>
const typename PetscSparseMatrix<TDataType>::Entry& PetscSparseMatrix<TDataType>::Entry::operator*= (const TDataType& value) const
{
    const PetscInt    petsc_i   = mRow;
    const PetscInt    petsc_j   = mCol;
    const PetscScalar new_value = C2P<TDataType> (static_cast<TDataType> (*this) * value);
    MatSetValues (( Mat ) mMatrix, 1, &petsc_i, 1, &petsc_j, &new_value, INSERT_VALUES);

    return *this;
}

template <typename TDataType>
const typename PetscSparseMatrix<TDataType>::Entry& PetscSparseMatrix<TDataType>::Entry::operator/= (const TDataType& value) const
{
    const PetscInt    petsc_i   = mRow;
    const PetscInt    petsc_j   = mCol;
    const PetscScalar new_value = static_cast<TDataType> (*this) / value;
    MatSetValues (( Mat ) mMatrix, 1, &petsc_i, 1, &petsc_j, &new_value, INSERT_VALUES);

    return *this;
}

template <typename TDataType>
PetscSparseMatrix<TDataType>::Entry::operator TDataType () const
{
    // first verify that the requested
    // element is actually locally
    // available
    const PetscInt petsc_i       = mRow;
    const PetscInt petsc_j       = mCol;
    const PetscInt petsc_i_begin = mMatrix.row_start ();
    const PetscInt petsc_i_end   = mMatrix.row_end ();
    const PetscInt petsc_j_begin = mMatrix.col_start ();
    const PetscInt petsc_j_end   = mMatrix.col_end ();

    TDataType value = 0.0;

    if (petsc_i >= petsc_i_begin && petsc_i < petsc_i_end)
    {
        if (petsc_j >= petsc_j_begin && petsc_j < petsc_j_end)
        {
            PetscScalar petsc_value;
            MatGetValues (( Mat ) mMatrix, 1, &petsc_i, 1, &petsc_j, &petsc_value);

            value = P2C<TDataType> (petsc_value);

            return value;
        }
    }

    // assert (0
    //         && "You are trying to access an element of a matrix "
    //            "that is neither a locally owned element nor a "
    //            "ghost element of the matrix.");
#ifdef SG_VERBOSE
    std::cerr << "You are trying to access an element of a matrix"
              << "that is neither a locally owned element nor a "
              << "ghost element of the matrix." << std::endl;
#endif

    value = get_matrix_value<TDataType> (( Mat ) mMatrix, petsc_i, petsc_j);

    return value;
}

//=============================================================================
// PetscSparseMatrix

template <typename TDataType>
PetscSparseMatrix<TDataType>::PetscSparseMatrix ()
{
    d = new PetscSparseMatrixPrivate ();
    // d->M = 0;
    // d->N = 0;
    d->mat = PETSC_NULLPTR;

    // CHKERRV (MatCreate (PETSC_COMM_WORLD, &d->mat));
    // CHKERRV (MatSetSizes (d->mat, d->m, d->n, d->M, d->N));
    // CHKERRV (MatSetUp (d->mat));
}

template <typename TDataType>
PetscSparseMatrix<TDataType>::PetscSparseMatrix (Sparsity* sp)
{
    d = new PetscSparseMatrixPrivate ();

    /*
     * Performance tuning note: For problems of substantial size,
     * preallocation of matrix memory is crucial for attaining good
     * performance. See the matrix chapter of the users manual for details.
     *
     * see https://petsc.org/release/manual/getting_started/
     */
    // todo: Perform the preallocation from the matrix shape.
    // ...
    CHKERRV (MatCreate (PETSC_COMM_WORLD, &d->mat));
    CHKERRV (MatSetSizes (d->mat, d->m, d->n, d->M, d->N));
    CHKERRV (MatSetUp (d->mat));
}

template <typename TDataType>
PetscSparseMatrix<TDataType>::PetscSparseMatrix (Index_t M, Index_t N)
{
    d = new PetscSparseMatrixPrivate ();
    d->M = M;
    d->N = N;

    CHKERRV (MatCreate (PETSC_COMM_WORLD, &d->mat));
    CHKERRV (MatSetSizes (d->mat, d->m, d->n, d->M, d->N));
    CHKERRV (MatSetUp (d->mat));
}

template <typename TDataType>
PetscSparseMatrix<TDataType>::PetscSparseMatrix (Mat A, bool inc_ref_count)
{
    d      = new PetscSparseMatrixPrivate ();
    d->mat = A;
    if (inc_ref_count)
    {
        PetscObjectReference (( PetscObject ) d->mat);
    }
}

template <typename TDataType>
PetscSparseMatrix<TDataType>::PetscSparseMatrix (const PetscSparseMatrix<TDataType>& A)
{
    d = new PetscSparseMatrixPrivate ();
    // Duplicates a matrix including the non-zero structure.
    CHKERRV (MatDuplicate (( Mat ) A, MAT_COPY_VALUES, &d->mat));
}

template <typename TDataType>
PetscSparseMatrix<TDataType>::PetscSparseMatrix (PetscSparseMatrix<TDataType>&& A)
{
    d      = new PetscSparseMatrixPrivate ();
    d->mat = ( Mat ) A;
    PetscObjectReference (( PetscObject ) d->mat);
}

template <typename TDataType>
PetscSparseMatrix<TDataType>::~PetscSparseMatrix ()
{
    if (d && d->mat)
    {
        // exceptionless_semiparallel_only ();

        // If we encounter an error here, print a warning but otherwise
        // keep going since we may be recovering from an exception.
        CHKERRV (MatDestroy (&d->mat));
        d->mat = nullptr;
    }

    if (d)
    {
        delete d;
    }
}

template <typename TDataType>
void PetscSparseMatrix<TDataType>::clear ()
{
    if (d && d->mat)
    {
        // exceptionless_semiparallel_only ();

        // If we encounter an error here, print a warning but otherwise
        // keep going since we may be recovering from an exception.
        CHKERRV (MatDestroy (&d->mat));
        d->mat = nullptr;
    }
}

template <typename TDataType>
void PetscSparseMatrix<TDataType>::close ()
{

}

template <typename TDataType>
PetscSparseMatrix<TDataType>& PetscSparseMatrix<TDataType>::operator= (const PetscSparseMatrix<TDataType>& A)
{
    // Duplicates a matrix including the non-zero structure.
    MatDuplicate (( Mat ) A, MAT_COPY_VALUES, &d->mat);
    return *this;
}

template <typename TDataType>
PetscSparseMatrix<TDataType>& PetscSparseMatrix<TDataType>::operator= (PetscSparseMatrix<TDataType>&& A)
{
    d->mat = ( Mat ) A;
    PetscObjectReference (( PetscObject ) d->mat);

    return *this;
}

template <typename TDataType>
PetscSparseMatrix<TDataType> PetscSparseMatrix<TDataType>::clone () const
{
    Mat mat;
    // MatDuplicate (d->mat, MAT_DO_NOT_COPY_VALUES, &mat);
    // MatCopy (d->mat, mat, SAME_NONZERO_PATTERN);
    MatDuplicate (d->mat, MAT_COPY_VALUES, &mat);
    PetscSparseMatrix<TDataType> Y (mat, true);
    MatDestroy (&mat);
    return Y;
}

template <typename TDataType>
Space<TDataType>* PetscSparseMatrix<TDataType>::getSpace ()
{
    static PetscSpace<TDataType> petsc_space;
    return &petsc_space;
}

template <typename TDataType>
PetscSparseMatrix<TDataType>::operator const Mat& () const
{
    return d->mat;
}

template <typename TDataType>
PetscSparseMatrix<TDataType>::operator Mat& ()
{
    return d->mat;
}

template <typename TDataType>
Mat PetscSparseMatrix<TDataType>::mat ()
{
    return d->mat;
}

template <typename TDataType>
Mat PetscSparseMatrix<TDataType>::mat () const
{
    return d->mat;
}

// template <typename TDataType>
// void PetscSparseMatrix<TDataType>::init (Sparsity* sp)
// {
//     CHKERRV (MatCreate (PETSC_COMM_WORLD, &d->mat));

//     /*
//      * The user could alternatively specify each processes’ number of local rows
//      * and columns using m and n.
//      */
//     CHKERRV (MatSetSizes (d->mat, d->m, d->n, d->M, d->N));

//     /*
//      * Performance tuning note: For problems of substantial size, preallocation of matrix
//      * memory is crucial for attaining good  performance. See the matrix chapter of the
//      * users manual for details.
//      */

//     this->mInitialized = true;
// }

// template <typename TDataType>
// void PetscSparseMatrix<TDataType>::clear ()
// {
//     if ((this->initialized ()) /*&& (this->_destroy_mat_on_exit)*/)
//     {
//         // exceptionless_semiparallel_only ();

//         // If we encounter an error here, print a warning but otherwise
//         // keep going since we may be recovering from an exception.
//         CHKERRV (MatDestroy (&d->mat));

//         // if (ierr)
//         //     libmesh_warning ("Warning: MatDestroy returned a non-zero error code which we ignored.");

//         this->mInitialized = false;
//     }
// }

template <typename TDataType>
void PetscSparseMatrix<TDataType>::beginAssembly ()
{
    CHKERRV (MatAssemblyBegin (d->mat, MAT_FINAL_ASSEMBLY));
}

template <typename TDataType>
void PetscSparseMatrix<TDataType>::endAssembly ()
{
    CHKERRV (MatAssemblyEnd (d->mat, MAT_FINAL_ASSEMBLY));
}

template <typename TDataType>
void PetscSparseMatrix<TDataType>::assemble ()
{
    CHKERRV (MatAssemblyBegin (d->mat, MAT_FINAL_ASSEMBLY));
    CHKERRV (MatAssemblyEnd (d->mat, MAT_FINAL_ASSEMBLY));
}

template <typename TDataType>
bool PetscSparseMatrix<TDataType>::assembled () const
{
    PetscBool has_assembled;
    CHKERRQ (MatAssembled (d->mat, &has_assembled));

    return has_assembled == PETSC_TRUE ? true : false;
}

template <typename TDataType>
void PetscSparseMatrix<TDataType>::resize (Index_t M, Index_t N)
{
    // has the exact dimensions?
    if (d->mat != PETSC_NULLPTR)
    {
        PetscInt M1, N1;
        CHKERRV (MatGetSize (d->mat, &M1, &N1));

        if (M1 == M && N1 == N)
        {
            return;
        }
    }

    if (d->mat != PETSC_NULLPTR)
    {
        CHKERRV (MatDestroy (&d->mat));
    }

    CHKERRV (MatCreate (PETSC_COMM_WORLD, &d->mat));
    CHKERRV (MatSetSizes (d->mat, PETSC_DECIDE, PETSC_DECIDE, M, N));
    CHKERRV (MatSetUp (d->mat));
    }

template <typename TDataType>
void PetscSparseMatrix<TDataType>::conservativeResize (Index_t M, Index_t N)
{
    // d->M = M;
    // d->N = N;

    if (d->mat == PETSC_NULLPTR)
    {
        CHKERRV (MatCreate (PETSC_COMM_WORLD, &d->mat));
        CHKERRV (MatSetSizes (d->mat, PETSC_DECIDE, PETSC_DECIDE, M, N));
        CHKERRV (MatSetUp (d->mat));

        return;
    }

    // once the vector or matrices sizes have been set and the matrices or vectors
    // are fully usable one cannot change the size of the matrices or vectors or
    // number of processors they live on. One may create new vectors and copy,
    // for example using VecScatterCreate(), the old values from the previous vector.
    Mat mat;
    CHKERRV (MatCreate (PETSC_COMM_WORLD, &mat));
    CHKERRV (MatSetSizes (mat, PETSC_DECIDE, PETSC_DECIDE, M, N));
    CHKERRV (MatSetUp (mat));

    CHKERRV (MatDestroy (&d->mat));
    d->mat = mat;
}

template <typename TDataType>
Index_t PetscSparseMatrix<TDataType>::rows () const
{
    PetscInt m, n;
    CHKERRQ (MatGetSize (d->mat, &m, &n));

    return m;
}

template <typename TDataType>
Index_t PetscSparseMatrix<TDataType>::cols () const
{
    PetscInt m, n;
    CHKERRQ (MatGetSize (d->mat, &m, &n));

    return n;
}

template <typename TDataType>
Index_t PetscSparseMatrix<TDataType>::row_start () const
{
    PetscInt start = 0, stop = 0;

    CHKERRQ (MatGetOwnershipRange (d->mat, &start, &stop));

    return static_cast<Index_t> (start);
}

template <typename TDataType>
Index_t PetscSparseMatrix<TDataType>::row_end () const
{
    PetscInt start = 0, stop = 0;

    CHKERRQ (MatGetOwnershipRange (d->mat, &start, &stop));

    return static_cast<Index_t> (stop);
}

template <typename TDataType>
Index_t PetscSparseMatrix<TDataType>::col_start () const
{
    PetscInt start = 0, stop = 0;

    CHKERRQ (MatGetOwnershipRangeColumn (d->mat, &start, &stop));

    return static_cast<Index_t> (start);
}

template <typename TDataType>
Index_t PetscSparseMatrix<TDataType>::col_end () const
{
    PetscInt start = 0, stop = 0;

    CHKERRQ (MatGetOwnershipRangeColumn (d->mat, &start, &stop));

    return static_cast<Index_t> (stop);
}

template <typename TDataType>
void PetscSparseMatrix<TDataType>::setZero ()
{
    CHKERRV (MatZeroEntries (d->mat));
}

template <typename TDataType>
void PetscSparseMatrix<TDataType>::setIdentity ()
{
    CHKERRV (MatZeroEntries (d->mat));

    PetscInt start = 0, stop = 0;

    CHKERRV (MatGetOwnershipRange (d->mat, &start, &stop));

    for (PetscInt i = start; i < stop; ++i)
    {
        PetscInt    petsc_i     = i;
        PetscScalar petsc_value = 1.0;
        CHKERRV (MatSetValues (d->mat, 1, &petsc_i, 1, &petsc_i, &petsc_value, INSERT_VALUES));
    }

    CHKERRV (MatAssemblyBegin (d->mat, MAT_FINAL_ASSEMBLY));
    CHKERRV (MatAssemblyEnd (d->mat, MAT_FINAL_ASSEMBLY));
}

template <typename TDataType>
void PetscSparseMatrix<TDataType>::set (const Index_t& i, const Index_t& j, const TDataType& value)
{
    PetscInt i_val = i, j_val = j;

    PetscScalar petsc_value = static_cast<PetscScalar> (value);
    CHKERRV (MatSetValues (d->mat, 1, &i_val, 1, &j_val, &petsc_value, INSERT_VALUES));
}

template <typename TDataType>
void PetscSparseMatrix<TDataType>::add (const Index_t& i, const Index_t& j, const TDataType& value)
{
    PetscInt i_val = i, j_val = j;

    PetscScalar petsc_value = static_cast<PetscScalar> (value);
    CHKERRV (MatSetValues (d->mat, 1, &i_val, 1, &j_val, &petsc_value, ADD_VALUES));
}

// template <typename TDataType>
// void PetscSparseMatrix<TDataType>::add (const DenseMatrix<TDataType>& dm, const IndexMap_t& rows, const IndexMap_t& cols)
// {}

// template <typename TDataType>
// void PetscSparseMatrix<TDataType>::add (const SparseMatrix<TDataType>& spm, const IndexMap_t& rows, const IndexMap_t& cols)
// {}

template <typename TDataType>
typename PetscSparseMatrix<TDataType>::Entry PetscSparseMatrix<TDataType>::operator() (const Index_t& i, const Index_t& j)
{
    Entry e (*this, i, j);
    return e;
}

template <typename TDataType>
TDataType PetscSparseMatrix<TDataType>::operator() (const Index_t& i, const Index_t& j) const
{
    // first verify that the requested
    // element is actually locally
    // available
    PetscInt petsc_i = i;
    PetscInt petsc_j = j;
    PetscInt petsc_i_begin, petsc_i_end;
    PetscInt petsc_j_begin, petsc_j_end;

    CHKERRQ (MatGetOwnershipRange (d->mat, &petsc_i_begin, &petsc_i_end));
    CHKERRQ (MatGetOwnershipRangeColumn (d->mat, &petsc_j_begin, &petsc_j_end));

    TDataType value = 0.0;

    if (petsc_i >= petsc_i_begin && petsc_i < petsc_i_end)
    {
        if (petsc_j >= petsc_j_begin && petsc_j < petsc_j_end)
        {
            PetscScalar petsc_value;
            CHKERRQ (MatGetValues (d->mat, 1, &petsc_i, 1, &petsc_j, &petsc_value));

            value = P2C<TDataType> (petsc_value);

            return value;
        }
    }

    // assert (0
    //         && "You are trying to access an element of a matrix "
    //            "that is neither a locally owned element nor a "
    //            "ghost element of the matrix.");
#ifdef SG_VERBOSE
    std::cerr << "You are trying to access an element of a matrix "
              << "that is neither a locally owned element nor a "
              << "ghost element of the matrix." << std::endl;
#endif

    value = get_matrix_value<TDataType> (d->mat, i, j);

    return value;
}

template <typename TDataType>
Index_t PetscSparseMatrix<TDataType>::nonZeros () const
{
    MatInfo info;
    CHKERRQ (MatGetInfo (d->mat, MAT_GLOBAL_SUM, &info));
    // PCHKERRQ (d->mat, ierr);
    return ( Index_t ) info.nz_used;
}

template <typename TDataType>
void PetscSparseMatrix<TDataType>::get_csr (std::vector<Index_t>&   indices_i,
                                            std::vector<Index_t>&   indices_j,
                                            std::vector<TDataType>& values) const
{
    PetscInt  shift     = 0;            // 0 or 1 indicating we want the indices starting at 0 or 1
    PetscBool symmetric = PETSC_FALSE;  // PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be symmetrized
    PetscBool inodecompressed =
        PETSC_FALSE;       // PETSC_TRUE or PETSC_FALSE indicating if the nonzero structure of the inodes or the nonzero elements is wanted.
    PetscInt        n;     // number of local rows in the(possibly compressed) matrix, use NULL if not needed
    const PetscInt* ia;    // the row pointers; that is ia[0] = 0, ia[row] = ia[row - 1] + number of elements in that row of the matrix, use
                           // NULL if not needed
    const PetscInt* ja;    // the column indices, use NULL if not needed
    PetscBool       done;  // indicates if the routine actually worked and returned appropriate ia[] and ja[] arrays;
                           // callers are responsible for handling the case when done == PETSC_FALSE and ia and ja are not set

    CHKERRV (MatGetRowIJ (d->mat, shift, symmetric, inodecompressed, &n, &ia, &ja, &done));

    indices_i.resize (n + 1);
    for (auto i = 0; i < n + 1; ++i)
    {
        indices_i[i] = ia[i];
    }

    Index_t nnz = indices_i[n];
    indices_j.resize (nnz);

    for (auto j = 0; j < nnz; ++j)
    {
        indices_j[j] = ja[j];
    }

    CHKERRV (MatRestoreRowIJ (d->mat, shift, symmetric, inodecompressed, &n, &ia, &ja, &done));

    // Copy the data
    values.resize (nnz);

    PetscInt start = 0, stop = 0;
    CHKERRV (MatGetOwnershipRange (d->mat, &start, &stop));

    Index_t nonzeros = 0;

    for (PetscInt i = start; i < stop; ++i)
    {
        const PetscScalar* petsc_row;

        PetscInt ncols = 0;

        CHKERRV (MatGetRow (d->mat, i, &ncols, PETSC_NULLPTR, &petsc_row));

        for (auto j = 0; j < ncols; ++j)
        {
            values[nonzeros++] = P2C<TDataType> (petsc_row[j]);
        }

        CHKERRV (MatRestoreRow (d->mat, i, &ncols, PETSC_NULLPTR, &petsc_row));
    }
}

template <typename TDataType>
void PetscSparseMatrix<TDataType>::print (std::ostream& out) const
{
    CHKERRV (MatView (d->mat, PETSC_VIEWER_STDOUT_WORLD));
}

template <typename TDataType>
void PetscSparseMatrix<TDataType>::create_submatrix (SparseMatrix<TDataType>& submatrix,
                                                     const GeneralIndexSet&   rows,
                                                     const GeneralIndexSet&   cols) const
{
    // Make sure the SparseMatrix passed in is really a PetscMatrix
    PetscSparseMatrix<TDataType>* petsc_submatrix = dynamic_cast<PetscSparseMatrix<TDataType>*> (&submatrix);

    // If we're not reusing submatrix and submatrix is already initialized
    // then we need to clear it, otherwise we get a memory leak.
    submatrix.clear ();

    MPI_Comm comm;
    CHKERRV (PetscObjectGetComm (( PetscObject ) d->mat, &comm));

    PetscInt petsc_submatrix_M = rows.size ();
    PetscInt petsc_submatrix_N = cols.size ();

    std::vector<PetscInt> petsc_submatrix_rows (petsc_submatrix_M);
    std::copy (rows.begin (), rows.end (), petsc_submatrix_rows.begin ());

    std::vector<PetscInt> petsc_submatrix_cols (petsc_submatrix_N);
    std::copy (cols.begin (), cols.end (), petsc_submatrix_cols.begin ());

    // Construct row and column index sets.
    IS isrow;
    CHKERRV (ISCreateGeneral (comm, petsc_submatrix_M, petsc_submatrix_rows.data (), PETSC_USE_POINTER, &isrow));

    IS iscol;
    CHKERRV (ISCreateGeneral (comm, petsc_submatrix_N, petsc_submatrix_cols.data (), PETSC_USE_POINTER, &iscol));

    // Extract submatrix
#if PETSC_VERSION_LT(3, 8, 0)
#   define CreateSubMatrix MatGetSubMatrix
#else
#   define CreateSubMatrix MatCreateSubMatrix
#endif

    // petsc_submatrix->resize (petsc_submatrix_M, petsc_submatrix_N);

    CHKERRV (CreateSubMatrix (d->mat, isrow, iscol, MAT_INITIAL_MATRIX, (&petsc_submatrix->d->mat)));

    CHKERRV (ISDestroy (&isrow));
    CHKERRV (ISDestroy (&iscol));
}

template <typename TDataType>
void PetscSparseMatrix<TDataType>::create_submatrix (SparseMatrix<TDataType>& submatrix,
                                                     const StrideIndexSet&    rows,
                                                     const StrideIndexSet&    cols) const
{
    // Make sure the SparseMatrix passed in is really a PetscMatrix
    PetscSparseMatrix<TDataType>* petsc_submatrix = dynamic_cast<PetscSparseMatrix<TDataType>*> (&submatrix);

    // If we're not reusing submatrix and submatrix is already initialized
    // then we need to clear it, otherwise we get a memory leak.
    submatrix.clear ();

    MPI_Comm comm;
    CHKERRV (PetscObjectGetComm (( PetscObject ) d->mat, &comm));

    PetscInt petsc_submatrix_M = rows.n;
    PetscInt petsc_submatrix_N = cols.n;

    // Construct row and column index sets.
    IS isrow;
    CHKERRV (ISCreateStride (comm, rows.n, rows.first, rows.step, &isrow));

    IS iscol;
    CHKERRV (ISCreateStride (comm, cols.n, cols.first, cols.step, &iscol));

    // Extract submatrix
#if PETSC_VERSION_LT(3, 8, 0)
#   define CreateSubMatrix MatGetSubMatrix
#else
#   define CreateSubMatrix MatCreateSubMatrix
#endif

    // petsc_submatrix->resize (petsc_submatrix_M, petsc_submatrix_N);

    // CHKERRV (CreateSubMatrix (d->mat, isrow, iscol, MAT_INITIAL_MATRIX, (&petsc_submatrix->d->mat)));
    CHKERRV (CreateSubMatrix (d->mat, isrow, PETSC_NULL, MAT_INITIAL_MATRIX, (&petsc_submatrix->d->mat)));

    CHKERRV (ISDestroy (&isrow));
    CHKERRV (ISDestroy (&iscol));
}

template <typename TDataType>
void PetscSparseMatrix<TDataType>::create_submatrix_nosort (SparseMatrix<TDataType>& submatrix,
                                                            const GeneralIndexSet&   rows,
                                                            const GeneralIndexSet&   cols) const
{
    // Make sure the SparseMatrix passed in is really a PetscMatrix
    PetscSparseMatrix<TDataType>* petsc_submatrix = dynamic_cast<PetscSparseMatrix<TDataType>*> (&submatrix);
    petsc_submatrix->resize (rows.size(), cols.size());
    CHKERRV (MatZeroEntries (petsc_submatrix->d->mat));

    MPI_Comm comm;
    PetscObjectGetComm (( PetscObject ) d->mat, &comm);

    PetscInt row_start = 0, row_stop = 0;
    CHKERRV (MatGetOwnershipRange (d->mat, &row_start, &row_stop));

    PetscInt           pc_ncols = 0;
    const PetscInt*    pc_cols;
    const PetscScalar* pc_vals;

    // data for creating the submatrix
    std::vector<PetscInt>    sub_cols;
    std::vector<PetscScalar> sub_vals;

    for (auto i = 0; i < rows.size (); ++i)
    {
        PetscInt sub_rid[] = { static_cast<PetscInt> (i) };
        PetscInt rid       = static_cast<PetscInt> (rows[i]);
        // only get value from local rows, and set values to the corresponding columns in the submatrix
        if (rows[i] >= row_start && rows[i] < row_stop)
        {
            // get one row of data from the original matrix
            CHKERRV (MatGetRow (d->mat, rid, &pc_ncols, &pc_cols, &pc_vals));
            // extract data from certain cols, save the indices and entries sub_cols and sub_vals
            for (auto j = 0; j < rows.size (); ++j)
            {
                for (unsigned int idx = 0; idx < static_cast<unsigned int> (pc_ncols); idx++)
                {
                    if (pc_cols[idx] == static_cast<PetscInt> (cols[j]))
                    {
                        sub_cols.push_back (static_cast<PetscInt> (j));
                        sub_vals.push_back (pc_vals[idx]);
                    }
                }
            }
            // set values
            CHKERRV (MatSetValues (petsc_submatrix->d->mat,
                                   1,
                                   sub_rid,
                                   static_cast<PetscInt> (sub_vals.size ()),
                                   sub_cols.data (),
                                   sub_vals.data (),
                                   INSERT_VALUES));
            CHKERRV (MatRestoreRow (d->mat, rid, &pc_ncols, &pc_cols, &pc_vals));
            // clear data for this row
            sub_cols.clear ();
            sub_vals.clear ();
        }
    }

    CHKERRV (MatAssemblyBegin (petsc_submatrix->d->mat, MAT_FINAL_ASSEMBLY));
    CHKERRV (MatAssemblyEnd (petsc_submatrix->d->mat, MAT_FINAL_ASSEMBLY));
}

template <typename TDataType>
void PetscSparseMatrix<TDataType>::create_submatrix_nosort (SparseMatrix<TDataType>& submatrix,
                                                            const StrideIndexSet&    stride_rows,
                                                            const StrideIndexSet&    stride_cols) const
{
    // the interested row/col index set
    GeneralIndexSet rows, cols;
    for (auto i = 0; i < stride_rows.n; ++i)
    {
        rows.push_back (stride_rows.first + i * stride_rows.step);
    }

    for (auto i = 0; i < stride_cols.n; ++i)
    {
        cols.push_back (stride_cols.first + i * stride_cols.step);
    }

    create_submatrix_nosort (submatrix, rows, cols);
}

template <typename TDataType>
void PetscSparseMatrix<TDataType>::get_diagonal (Vector<TDataType>& dest) const
{
    // Make sure the NumericVector passed in is really a PetscVector
    PetscVector<TDataType>* petsc_dest = dynamic_cast<PetscVector<TDataType>*> (&dest);

    // Needs a const_cast since PETSc does not work with const.
    CHKERRV (MatGetDiagonal (d->mat, petsc_dest->vec ()));
}

template <typename TDataType>
void PetscSparseMatrix<TDataType>::get_transpose (SparseMatrix<TDataType>& dest) const
{
    // Make sure the SparseMatrix passed in is really a PetscMatrix
    PetscSparseMatrix<TDataType>* petsc_dest = dynamic_cast<PetscSparseMatrix<TDataType>*> (&dest);

    // If we aren't reusing the matrix then need to clear dest,
    // otherwise we get a memory leak
    if (petsc_dest != this)
        dest.clear ();

    if (petsc_dest == this)
    // The MAT_REUSE_MATRIX flag was replaced by MAT_INPLACE_MATRIX
    // in PETSc 3.7.0
#if PETSC_VERSION_LT(3, 7, 0)
        CHKERRV (MatTranspose (d->mat, MAT_REUSE_MATRIX, &petsc_dest->d->mat));
#else
        CHKERRV (MatTranspose (d->mat, MAT_INPLACE_MATRIX, &petsc_dest->d->mat));
#endif
    else
        CHKERRV (MatTranspose (d->mat, MAT_INITIAL_MATRIX, &petsc_dest->d->mat));
}

template <typename TDataType>
void PetscSparseMatrix<TDataType>::get_inverse (SparseMatrix<TDataType>& dest) const
{
    // Make sure the SparseMatrix passed in is really a PetscMatrix
    PetscSparseMatrix<TDataType>* petsc_dest = dynamic_cast<PetscSparseMatrix<TDataType>*> (&dest);
    MPI_Comm                      comm;
    PetscObjectGetComm (( PetscObject ) d->mat, &comm);

    PetscInt M, N;
    CHKERRV (MatGetSize (d->mat, &M, &N));

    PetscInt m, n;
    CHKERRV (MatGetLocalSize (d->mat, &m, &n));

    PetscInt row_start = 0, row_stop = 0;
    CHKERRV (MatGetOwnershipRange (d->mat, &row_start, &row_stop));

    /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                Create the linear solver and set various options
    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    KSP ksp; /* linear solver context */
    PC  pc;  /* preconditioner context */
    CHKERRV (KSPCreate (comm, &ksp));

    // KSPSetUp (ksp);

    /*
     Set runtime options, e.g.,
         -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
     These options will override those specified above as long as
     KSPSetFromOptions() is called _after_ any other customization
     routines.
     */
    CHKERRV (KSPSetFromOptions (ksp));

    /*
     Set operators. Here the matrix that defines the linear system
     also serves as the matrix that defines the preconditioner.
    */
    CHKERRV (KSPSetOperators (ksp, d->mat, d->mat));

    /*
     Set linear solver defaults for this problem (optional).
     - By extracting the KSP and PC contexts from the KSP context,
       we can then directly call any KSP and PC routines to set
       various options.
     - The following four statements are optional; all of these
       parameters could alternatively be specified at runtime via
       KSPSetFromOptions();
    */
    CHKERRV (KSPGetPC (ksp, &pc));
    CHKERRV (PCSetType (pc, PCJACOBI));
    CHKERRV (KSPSetTolerances (ksp, 1.e-5, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));

    // generate the identity matrix...
    // see https://petsc.org/release/manualpages/Mat/MatCreateDense/
    Mat E;
    CHKERRV (MatCreateDense (comm, m, n, M, N, NULL, &E));
    CHKERRV (MatZeroEntries (E));
    for (auto i = row_start; i < row_stop; ++i)
    {
        CHKERRV (MatSetValue (E, i, i, 1.0, INSERT_VALUES));
    }
    CHKERRV (MatAssemblyBegin (E, MAT_FINAL_ASSEMBLY));
    CHKERRV (MatAssemblyEnd (E, MAT_FINAL_ASSEMBLY));

    // the inverse matrix usually dense matrix
    Mat X;
    CHKERRV (MatDuplicate (E, MAT_DO_NOT_COPY_VALUES, &X));

    // Solves a linear system with multiple right-hand sides stored as a MATDENSE.
    // see https://petsc.org/release/manualpages/KSP/KSPMatSolve/
    CHKERRV (KSPMatSolve (ksp, E, X));

    // clear the petsc resource allocated
    petsc_dest->clear ();
    // update the petec handle
    petsc_dest->d->mat = X;

    /*
    Free work space.  All PETSc objects should be destroyed when they
    are no longer needed.
    */
    CHKERRV (KSPDestroy (&ksp));
    CHKERRV (MatDestroy (&E));
}

template <typename TDataType>
void PetscSparseMatrix<TDataType>::scale (const TDataType s)
{
    PetscScalar a = C2P<TDataType> (s);
    CHKERRV (MatScale (d->mat, a));
}

template <typename TDataType>
void PetscSparseMatrix<TDataType>::matrix_matrix_mult (SparseMatrix<TDataType>& X_in, SparseMatrix<TDataType>& Y_out)
{
    /*
     * Set reuse = true if this->_mat and X have the same nonzero pattern as before,
     * and Y is obtained from a previous call to this function with reuse = false
     */
    bool                          reuse = false;

    PetscSparseMatrix<TDataType>* X     = dynamic_cast<PetscSparseMatrix<TDataType>*> (&X_in);
    PetscSparseMatrix<TDataType>* Y = dynamic_cast<PetscSparseMatrix<TDataType>*> (&Y_out);

    if (reuse)
        CHKERRV (MatMatMult (d->mat, X->d->mat, MAT_REUSE_MATRIX, PETSC_DEFAULT, &Y->d->mat));
    else
    {
        Y->clear ();
        CHKERRV (MatMatMult (d->mat, X->d->mat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Y->d->mat));
    }
}

//=============================================================================
// Maintain compatibilities with Eigen

// template <typename TDataType>
// void PetscVector<TDataType>::add (const Eigen::VectorX<TDataType>& vec, const IndexMap_t& indices)
// {
//     std::cerr << "PetscVector<TDataType>::add NOT IMPLEMENTED" << std::endl;
// }

// template <typename TDataType>
// void PetscVector<TDataType>::get (Eigen::VectorX<TDataType>& vec, const IndexMap_t& indices)
// {
//     std::cerr << "PetscVector<TDataType>::get NOT IMPLEMENTED" << std::endl;
// }

// template <typename TDataType>
// TDataType* PetscVector<TDataType>::data ()
// {
//     // assert (0 && "PetscVector<TDataType>::data NOT IMPLEMENTED");
//     std::cerr << "PetscVector<TDataType>::data NOT IMPLEMENTED" << std::endl;
//     return ( TDataType* ) nullptr;
// }

// template <typename TDataType>
// const TDataType* PetscVector<TDataType>::data () const
// {
//     // assert (0 && "PetscVector<TDataType>::data NOT IMPLEMENTED");
//     std::cerr << "PetscVector<TDataType>::data NOT IMPLEMENTED" << std::endl;
//     return ( TDataType* ) nullptr;
// }

// template <typename TDataType>
// TDataType PetscVector<TDataType>::dot (const PetscVector<TDataType>& other) const
// {
//     PetscScalar petsc_value = 0.0;
//     CHKERRQ (VecDot (d->vec, other.d->vec, &petsc_value));
//     auto value = P2C<TDataType> (petsc_value);
//     return value;
// }

// template <typename TDataType>
// void PetscSparseMatrix<TDataType>::add (const Eigen::SparseMatrix<TDataType>& spm, const IndexMap_t& rows, const IndexMap_t& cols)
// {
//     std::cerr << "PetscSparseMatrix<TDataType>::add NOT IMPLEMENTED" << std::endl;
// }

// template <typename TDataType>
// void PetscSparseMatrix<TDataType>::resize (Index_t M, Index_t N)
// {
//     d->M = M;
//     d->N = N;


//     CHKERRV (MatSetSizes (d->mat, d->m, d->M, d->n, d->N));
// }

// template <typename TDataType>
// Index_t PetscSparseMatrix<TDataType>::nonZeros () const
// {
//     MatInfo info;
//     CHKERRQ (MatGetInfo (d->mat, MAT_GLOBAL_SUM, &info));
//     // PCHKERRQ (d->mat, ierr);
//     return ( Index_t ) info.nz_used;
// }

// template <typename TDataType>
// const Index_t* PetscSparseMatrix<TDataType>::outerIndexPtr () const
// {
//     // assert (0 && "PetscMatrix<TDataType>::outerIndexPtr NOT IMPLEMENTED");
//     std::cerr << "PetscMatrix<TDataType>::outerIndexPtr NOT IMPLEMENTED" << std::endl;
//     return ( Index_t* ) nullptr;
// }

// template <typename TDataType>
// Index_t* PetscSparseMatrix<TDataType>::innerIndexPtr ()
// {
//     // assert (0 && "PetscMatrix<TDataType>::innerIndexPtr NOT IMPLEMENTED");
//     std::cerr << "PetscMatrix<TDataType>::innerIndexPtr NOT IMPLEMENTED" << std::endl;
//     return ( Index_t* ) nullptr;
// }

// template <typename TDataType>
// const Index_t* PetscSparseMatrix<TDataType>::innerIndexPtr () const
// {
//     // assert (0 && "PetscMatrix<TDataType>::innerIndexPtr NOT IMPLEMENTED");
//     std::cerr << "PetscMatrix<TDataType>::innerIndexPtr NOT IMPLEMENTED" << std::endl;
//     return ( Index_t* ) nullptr;
// }

// template <typename TDataType>
// TDataType* PetscSparseMatrix<TDataType>::valuePtr ()
// {
//     // assert (0 && "PetscMatrix<TDataType>::valuePtr NOT IMPLEMENTED");
//     std::cerr << "PetscSparseMatrix<TDataType>::valuePtr NOT IMPLEMENTED" << std::endl;
//     return ( TDataType* ) nullptr;
// }

// template <typename TDataType>
// const TDataType* PetscSparseMatrix<TDataType>::valuePtr () const
// {
//     // assert (0 && "PetscMatrix<TDataType>::valuePtr NOT IMPLEMENTED");
//     std::cerr << "PetscSparseMatrix<TDataType>::valuePtr NOT IMPLEMENTED" << std::endl;
//     return ( TDataType* ) nullptr;
// }

// template <typename TDataType>
// void PetscSparseMatrix<TDataType>::makeCompressed ()
// {
//     assemble ();
// }

//=============================================================================
// Maintain compatibilities with Eigen

// template <typename TDataType>
// void PetscVector<TDataType>::add (const Eigen::VectorX<TDataType>& vec, const IndexMap_t& indices)
// {
//     std::cerr << "PetscVector<TDataType>::add NOT IMPLEMENTED" << std::endl;
// }

// template <typename TDataType>
// void PetscVector<TDataType>::get (Eigen::VectorX<TDataType>& vec, const IndexMap_t& indices)
// {
//     std::cerr << "PetscVector<TDataType>::get NOT IMPLEMENTED" << std::endl;
// }

// template <typename TDataType>
// TDataType* PetscVector<TDataType>::data ()
// {
//     // assert (0 && "PetscVector<TDataType>::data NOT IMPLEMENTED");
//     std::cerr << "PetscVector<TDataType>::data NOT IMPLEMENTED" << std::endl;
//     return ( TDataType* ) nullptr;
// }

// template <typename TDataType>
// const TDataType* PetscVector<TDataType>::data () const
// {
//     // assert (0 && "PetscVector<TDataType>::data NOT IMPLEMENTED");
//     std::cerr << "PetscVector<TDataType>::data NOT IMPLEMENTED" << std::endl;
//     return ( TDataType* ) nullptr;
// }

// template <typename TDataType>
// TDataType PetscVector<TDataType>::dot (const PetscVector<TDataType>& other) const
// {
//     PetscScalar petsc_value = 0.0;
//     CHKERRQ (VecDot (d->vec, other.d->vec, &petsc_value));
//     auto value = P2C<TDataType> (petsc_value);
//     return value;
// }

// template <typename TDataType>
// void PetscSparseMatrix<TDataType>::add (const Eigen::SparseMatrix<TDataType>& spm, const IndexMap_t& rows, const IndexMap_t& cols)
// {
//     std::cerr << "PetscSparseMatrix<TDataType>::add NOT IMPLEMENTED" << std::endl;
// }

// template <typename TDataType>
// void PetscSparseMatrix<TDataType>::resize (Index_t M, Index_t N)
// {
//     d->M = M;
//     d->N = N;


//     CHKERRV (MatSetSizes (d->mat, d->m, d->M, d->n, d->N));
// }

// template <typename TDataType>
// Index_t PetscSparseMatrix<TDataType>::nonZeros () const
// {
//     MatInfo info;
//     CHKERRQ (MatGetInfo (d->mat, MAT_GLOBAL_SUM, &info));
//     // PCHKERRQ (d->mat, ierr);
//     return ( Index_t ) info.nz_used;
// }

// template <typename TDataType>
// const Index_t* PetscSparseMatrix<TDataType>::outerIndexPtr () const
// {
//     // assert (0 && "PetscMatrix<TDataType>::outerIndexPtr NOT IMPLEMENTED");
//     std::cerr << "PetscMatrix<TDataType>::outerIndexPtr NOT IMPLEMENTED" << std::endl;
//     return ( Index_t* ) nullptr;
// }

// template <typename TDataType>
// Index_t* PetscSparseMatrix<TDataType>::innerIndexPtr ()
// {
//     // assert (0 && "PetscMatrix<TDataType>::innerIndexPtr NOT IMPLEMENTED");
//     std::cerr << "PetscMatrix<TDataType>::innerIndexPtr NOT IMPLEMENTED" << std::endl;
//     return ( Index_t* ) nullptr;
// }

// template <typename TDataType>
// const Index_t* PetscSparseMatrix<TDataType>::innerIndexPtr () const
// {
//     // assert (0 && "PetscMatrix<TDataType>::innerIndexPtr NOT IMPLEMENTED");
//     std::cerr << "PetscMatrix<TDataType>::innerIndexPtr NOT IMPLEMENTED" << std::endl;
//     return ( Index_t* ) nullptr;
// }

// template <typename TDataType>
// TDataType* PetscSparseMatrix<TDataType>::valuePtr ()
// {
//     // assert (0 && "PetscMatrix<TDataType>::valuePtr NOT IMPLEMENTED");
//     std::cerr << "PetscSparseMatrix<TDataType>::valuePtr NOT IMPLEMENTED" << std::endl;
//     return ( TDataType* ) nullptr;
// }

// template <typename TDataType>
// const TDataType* PetscSparseMatrix<TDataType>::valuePtr () const
// {
//     // assert (0 && "PetscMatrix<TDataType>::valuePtr NOT IMPLEMENTED");
//     std::cerr << "PetscSparseMatrix<TDataType>::valuePtr NOT IMPLEMENTED" << std::endl;
//     return ( TDataType* ) nullptr;
// }

// template <typename TDataType>
// void PetscSparseMatrix<TDataType>::makeCompressed ()
// {
//     assemble ();
// }

//=============================================================================
// Operators

SG_NAMESPACE_OPEN
namespace Algebra
{
    template <typename TDataType>
    PetscVector<TDataType> operator+ (const PetscVector<TDataType>& x, const PetscVector<TDataType>& y)
    {
        PetscScalar a = C2P<TDataType> (1.0);
        auto        w = x.clone ();
        VecWAXPY (( Vec ) w, a, ( Vec ) x, ( Vec ) y);

        return w;
    }

    template <typename TDataType>
    PetscVector<TDataType> operator- (const PetscVector<TDataType>& x, const PetscVector<TDataType>& y)
    {
        PetscScalar a = C2P<TDataType> (-1.0);
        auto        w = x.clone ();
        VecWAXPY (( Vec ) w, a, ( Vec ) y, ( Vec ) x);

        return w;
    }

    template <typename TDataType>
    PetscVector<TDataType> operator* (const TDataType& s, const PetscVector<TDataType>& x)
    {
        PetscScalar a = C2P<TDataType> (s);
        auto        y = x.clone ();
        VecScale (( Vec ) y, a);

        return y;
    }

    template <typename TDataType>
    TDataType operator* (const PetscVector<TDataType>& x, const PetscVector<TDataType>& y)
    {
        PetscScalar petsc_value = 0.0;
        CHKERRQ (VecDot (( Vec ) x, ( Vec ) y, &petsc_value));
        auto value = P2C<TDataType> (petsc_value);

        return value;
    }

    template <typename TDataType>
    PetscVector<TDataType> operator/ (const PetscVector<TDataType>& x, const TDataType& s)
    {
        PetscScalar            a = C2P<TDataType> (1.0 / s);
        PetscVector<TDataType> y = x.clone ();
        VecScale (( Vec ) y, a);

        return y;
    }

    template <typename TDataType>
    PetscSparseMatrix<TDataType> operator+ (const PetscSparseMatrix<TDataType>& X, const PetscSparseMatrix<TDataType>& Y)
    {
        PetscScalar a = C2P<TDataType> (1.0);
        auto        W = Y.clone ();
        MatAYPX (( Mat ) W, a, ( Mat ) X, SAME_NONZERO_PATTERN);

        return W;
    }

    template <typename TDataType>
    PetscSparseMatrix<TDataType> operator- (const PetscSparseMatrix<TDataType>& X, const PetscSparseMatrix<TDataType>& Y)
    {
        PetscScalar a = C2P<TDataType> (-1.0);
        auto        W = Y.clone ();
        MatAYPX (( Mat ) W, a, ( Mat ) X, SAME_NONZERO_PATTERN);

        return W;
    }

    template <typename TDataType>
    PetscSparseMatrix<TDataType> operator* (const TDataType& s, const PetscSparseMatrix<TDataType>& A)
    {
        PetscScalar a = C2P<TDataType> (s);
        auto        W = A.clone ();
        MatScale (( Mat ) W, a);

        return W;
    }

    template <typename TDataType>
    PetscVector<TDataType> operator* (const PetscSparseMatrix<TDataType>& A, const PetscVector<TDataType>& x)
    {
        auto y = x.clone ();
        MatMult (( Mat ) A, ( Vec ) x, ( Vec ) y);

        return y;
    }

    template <typename TDataType>
    TDataType quadratic (const PetscSparseMatrix<TDataType>& A, const PetscVector<TDataType>& x)
    {
        auto y = x.clone ();
        MatMult (( Mat ) A, ( Vec ) x, ( Vec ) y);

        PetscScalar petsc_value = 0.0;
        CHKERRQ (VecDot (( Vec ) x, ( Vec ) y, &petsc_value));
        auto value = P2C<TDataType> (petsc_value);

        return value;
    }

    template <typename TDataType>
    std::ostream& operator<< (std::ostream& out, const PetscVector<TDataType>& v)
    {
        v.print (out);
        return out;
    }

    template <typename TDataType>
    std::ostream& operator<< (std::ostream& out, const PetscDenseMatrix<TDataType>& m)
    {
        m.print (out);
        return out;
    }

    template <typename TDataType>
    std::ostream& operator<< (std::ostream& out, const PetscSparseMatrix<TDataType>& m)
    {
        m.print (out);
        return out;
    }
}  // namespace Algebra
SG_NAMESPACE_CLOSE

//=============================================================================
// Explicit instantiations

SG_NAMESPACE_OPEN
namespace Algebra
{
    //------------------------------------------------------------------
    // Real_t
    template class ALGEBRA_EXPORT PetscVector<Real_t>;
    template class ALGEBRA_EXPORT PetscDenseMatrix<Real_t>;
    template class ALGEBRA_EXPORT PetscSparseMatrix<Real_t>;
    template class ALGEBRA_EXPORT PetscSpace<Real_t>;

    template ALGEBRA_EXPORT PetscVector<Real_t> operator+ <Real_t> (const PetscVector<Real_t>& x, const PetscVector<Real_t>& y);
    template ALGEBRA_EXPORT PetscVector<Real_t> operator- <Real_t> (const PetscVector<Real_t>&, const PetscVector<Real_t>&);
    template ALGEBRA_EXPORT PetscVector<Real_t> operator* <Real_t> (const Real_t& s, const PetscVector<Real_t>& x);
    template ALGEBRA_EXPORT Real_t operator* <Real_t> (const PetscVector<Real_t>& x, const PetscVector<Real_t>& y);
    template ALGEBRA_EXPORT PetscVector<Real_t> operator/ <Real_t> (const PetscVector<Real_t>& x, const Real_t& s);

    template ALGEBRA_EXPORT PetscSparseMatrix<Real_t> operator+ <Real_t> (const PetscSparseMatrix<Real_t>&,
                                                                          const PetscSparseMatrix<Real_t>&);
    template ALGEBRA_EXPORT PetscSparseMatrix<Real_t> operator- <Real_t> (const PetscSparseMatrix<Real_t>&,
                                                                          const PetscSparseMatrix<Real_t>&);
    template ALGEBRA_EXPORT PetscSparseMatrix<Real_t> operator* <Real_t> (const Real_t&, const PetscSparseMatrix<Real_t>&);
    template ALGEBRA_EXPORT PetscVector<Real_t> operator* <Real_t> (const PetscSparseMatrix<Real_t>&, const PetscVector<Real_t>&);

    template ALGEBRA_EXPORT Real_t quadratic<Real_t> (const PetscSparseMatrix<Real_t>&, const PetscVector<Real_t>&);

    template ALGEBRA_EXPORT std::ostream& operator<< <Real_t> (std::ostream& s, const PetscVector<Real_t>& v);
    template ALGEBRA_EXPORT std::ostream& operator<< <Real_t> (std::ostream& s, const PetscDenseMatrix<Real_t>& m);
    template ALGEBRA_EXPORT std::ostream& operator<< <Real_t> (std::ostream& s, const PetscSparseMatrix<Real_t>& m);

    //------------------------------------------------------------------
    // Complex_t
    template class ALGEBRA_EXPORT PetscVector<Complex_t>;
    template class ALGEBRA_EXPORT PetscDenseMatrix<Complex_t>;
    template class ALGEBRA_EXPORT PetscSparseMatrix<Complex_t>;
    template class ALGEBRA_EXPORT PetscSpace<Complex_t>;

    template ALGEBRA_EXPORT PetscVector<Complex_t> operator+ <Complex_t> (const PetscVector<Complex_t>& x, const PetscVector<Complex_t>& y);
    template ALGEBRA_EXPORT PetscVector<Complex_t> operator- <Complex_t> (const PetscVector<Complex_t>&, const PetscVector<Complex_t>&);
    template ALGEBRA_EXPORT PetscVector<Complex_t> operator* <Complex_t> (const Complex_t& s, const PetscVector<Complex_t>& x);
    template ALGEBRA_EXPORT Complex_t operator* <Complex_t> (const PetscVector<Complex_t>& x, const PetscVector<Complex_t>& y);
    template ALGEBRA_EXPORT PetscVector<Complex_t> operator/ <Complex_t> (const PetscVector<Complex_t>& x, const Complex_t& s);

    template ALGEBRA_EXPORT PetscSparseMatrix<Complex_t> operator+ <Complex_t> (const PetscSparseMatrix<Complex_t>&,
                                                                                const PetscSparseMatrix<Complex_t>&);
    template ALGEBRA_EXPORT PetscSparseMatrix<Complex_t> operator- <Complex_t> (const PetscSparseMatrix<Complex_t>&,
                                                                                const PetscSparseMatrix<Complex_t>&);
    template ALGEBRA_EXPORT PetscSparseMatrix<Complex_t> operator* <Complex_t> (const Complex_t&, const PetscSparseMatrix<Complex_t>&);
    template ALGEBRA_EXPORT PetscVector<Complex_t> operator* <Complex_t> (const PetscSparseMatrix<Complex_t>&,
                                                                          const PetscVector<Complex_t>&);

    template ALGEBRA_EXPORT Complex_t quadratic<Complex_t> (const PetscSparseMatrix<Complex_t>&, const PetscVector<Complex_t>&);

    template ALGEBRA_EXPORT std::ostream& operator<< <Complex_t> (std::ostream& s, const PetscVector<Complex_t>& m);
    template ALGEBRA_EXPORT std::ostream& operator<< <Complex_t> (std::ostream& s, const PetscDenseMatrix<Complex_t>& m);
    template ALGEBRA_EXPORT std::ostream& operator<< <Complex_t> (std::ostream& s, const PetscSparseMatrix<Complex_t>& m);

}  // namespace Algebra
SG_NAMESPACE_CLOSE

//=============================================================================
// Explicit specilizations

SG_NAMESPACE_OPEN
namespace Algebra
{
    //------------------------------------------------------------------
    // PetscVector<TDataType>

    //------------------------------------------------------------------
    // PetscSparseMatrix<TDataType>

}  // namespace Algebra
SG_NAMESPACE_CLOSE